
/*
 * main.c
 *
 *  Created on: Mar 29, 2010
 *      Author: nguyencao
 */

#include "xcun.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int m, col, b, n, pm, start, dcta, won, tdlim, pplim, del, LMlim, LMmod, LMres, LM, limit;
int tbl_first; /* start n of the whole table */
int dta[50];
INTEGER ppart[1000];
INTEGER ppartL[1000];
INTEGER ppartM[1000];
FILE* file;

int data[100][30];		// store lines read from xcun.d
int ntimes;

int main(int argc, char* argv[]) {

	int i, j, dctb, dtb[50],t,count1,count2;
	int islm;

	initialize();
	for(t=0;t<ntimes;t++) {
		initialize2(t);
		/*
		 * Begin main loop
		 * */
		for(; n<= LMlim; n+=del){
			/* The indicator of having LM's */
			islm = (LM ==1) && ((n%LMmod)==LMres);
			if(n>limit && !islm)
				continue;
			/* form a list of divisors of n whose codivisors are odd */
			/* dcta count the forward list, dctb count the backward list */
			dcta = 0;
			dctb = 0;

			if( (n%2 == 1)  && (n > 1) && !(b==2 && pm<0)){
				dcta = 1;
				dta[0]=1;
			}
			/* i =  divisor of n*/
			/* j = codivisor of i*/
			for(i=2; i<=n; i++){
				j = n/i;
				if(j < i) break;
				if( n != i*j) continue; /* j is not a divisor of n */

				if((i%2) == 1){
					dtb[dctb] = j;
					dctb++;
				}
				if(j <= i) break;
				if((j%2) == 1){
					dta[dcta] = i;
					dcta++;
				}
			}
			/* Copy the reverse list (odd co-divisors) onto the forward list */
			if(dctb > 0){
				for(i=0; i<dctb; i++){
					dta[dcta] = dtb[dctb-1-i];
					dcta++;
				}
			}
			/* Write n to stdout (main), 2 blanks after n */
			putdec2(STDOUT, n, 5);
			if(b != 2) {
				putc(' ', STDOUT);
				putc(' ', STDOUT);
				col = 6;
			} else {
				if(islm) {
					col = 6;
				} else {
					putc(' ', STDOUT);
					putc(' ', STDOUT);
					col = 6;
				}
			}

			if(islm){/* If n has LM's */
				dctb =0;
				for(i=0; i<dcta; i++){
					j = jacobi(b, n/dta[i]);
					if(j == 0){
						dctb++;
						if(b != 2) {
							if(dctb == 1){
								putc('(', STDOUT);
							}
							else{
								putc(',', STDOUT);
							}
						}

						if(b!=2) {
							if( (dta[i] % LMmod) == LMres){
								if(dta[i] != 3 && dta[i] != 2) {
									putdec2(STDOUT, dta[i], 1);
									putc('L', STDOUT);
									putc(',', STDOUT);
									putdec2(STDOUT, dta[i], 1);
									putc('M', STDOUT);
								} else {
									putdec2(STDOUT, dta[i], 1);
								}
							} else {
								putdec2(STDOUT, dta[i], 1);
							}
						}
					}
					/* Either set sign or kill it */
					dta[i] = j*dta[i];
				}
				if(dctb > 0 && b != 2) {
					putc(')', STDOUT);
					putc(' ', STDOUT);
				}

				/* Write L.M */
				if(b != 2) {
					putc('L', STDOUT);
					putc('.', STDOUT);
					putc('M', STDOUT);
					putc(NEWLINE, STDOUT);

					for(i=0; i<5; i++){
						putc(' ', STDOUT);
					}
				}
				putc('L', STDOUT);
				putc(' ', STDOUT);
				/* Count non-zero entries in dta */
				count1 = 0;
				for(i=0;i<dcta;i++) {
					if(dta[i] != 0)
						count1++;
				}
				j=0;
				count2 = 0;
				for(i=0; i<dcta; i++){
					if(dta[i] != 0){
						count2++;
						j++;
						if(abs(dta[i]) != 3 && abs(dta[i]) != 2)  {
							if(j == 1) {
								putc('(', STDOUT);
								col++;
							}
							col = col+idec(abs(dta[i]));
							if(dta[i]>0) {
								putc('L', STDOUT);
								col++;
							} else {
								putc('M', STDOUT);
								col++;
							}
							if(count2<count1){
								putc(',', STDOUT);
								col++;
							}
							else if(count1 == count2) {
								putc(')', STDOUT);
								col++;
							}
						} else {
							if(dta[i] < 0) {
								if(j == 1) {
									putc('(', STDOUT);
									col++;
								}
								col = col+idec(abs(dta[i]));
								if(count2<count1){
									putc(',', STDOUT);
									col++;
								}
								else if(count1 == count2){
									putc(')', STDOUT);
									col++;
								}
							} else {
								// 3L or 2L, do nothing
								if(count1 == 1) {
									// do nothing
								} else {
									if(j == 1) {
										putc('(', STDOUT);
										col++;
									}
									if(count1 == count2) {
										putc(')', STDOUT);
										col++;
									}
								}
							}
						}
					}
				}
				putc(' ', STDOUT);
				col++;
				proc(1);	/* L */

				for(i=0; i<5; i++){
					putc(' ', STDOUT);
				}

				putc('M', STDOUT);
				putc(' ', STDOUT);
				col=6;
				/* Count non-zero entries in dta */
				j=0;
				count2 = 0;
				for(i=0; i<dcta; i++){
					if(dta[i] != 0){
						count2++;
						j++;

						if(abs(dta[i]) != 3 && abs(dta[i]) != 2)  {
							if(j == 1) {
								putc('(', STDOUT);
								col++;
							}
							col = col+idec(abs(dta[i]));
							if(dta[i]<0) {
								putc('L', STDOUT);
								col++;
							} else {
								putc('M', STDOUT);
								col++;
							}
							if(count2<count1){
								putc(',', STDOUT);
								col++;
							}
							else if(count1 == count2) {
								putc(')', STDOUT);
								col++;
							}
						} else {
							if(dta[i] > 0) {
								if(j == 1) {
									putc('(', STDOUT);
									col++;
								}
								col = col+idec(abs(dta[i]));
								if(count2<count1){
									putc(',', STDOUT);
									col++;
								}
								else if(count1 == count2){
									putc(')', STDOUT);
									col++;
								}
							} else {
								// 3L or 2L, do nothing
								if(count1 == 1) {
									// do nothing
								} else {
									if(j == 1) {
										putc('(', STDOUT);
										col++;
									}
									if(count1 == count2) {
										putc(')', STDOUT);
										col++;
									}
								}
							}
						}
					}
				}
				putc(' ', STDOUT);
				col++;
				proc(-1);	/* M */
			} else{/* Else n does not have LM's, just dump the ( ) */
				if(dcta > 0){
					putc('(', STDOUT);
					col++;
					for(i=0; i<dcta; i++){
						if(i>=1){
							putc(',', STDOUT);
							col++;
						}
						/* col = col + length of divisor dta[i]+ space of '(' */
						col = col + idec(dta[i]);
					}
					putc(')', STDOUT);
					col++; 	/* space of ')' */
				}
				putc(' ', STDOUT);
				col++;
				proc(0);
			}
		}

		/*  end main loop
		 */
	}
	return 0;
}

void proc(int m) {
	int t, i, j, k;
	int pp_idx; /* index to access ppart[] */
	INTEGER	 ppt, tmp, quot, rem;
	char *bigdig;
	char *err_msg;
	init(ppt);
	init(tmp);
	init(quot);
	init(rem);
	if(m == 0){/* Not LM part */
		mpow(b, n, ppt);
		add2(ppt,won,ppt);
		if(dcta > 0){	/* Divide ppt by divisors in (), to get the part after ) */
			for(i=0; i< dcta; i++){
				pp_idx = (dta[i]- tbl_first)/del;
				mcopy(tmp, ppart[pp_idx]);	/* tmp = ppart[pp_idx] */
				div3(ppt, tmp, quot, rem);
				if(compare2(rem, 0) != 0){
					err_msg = "error dividing cofactor";
					error(err_msg);
				}
				mcopy(ppt, quot);
			}
		}
		if(n <= pplim){
			pp_idx = (n- tbl_first)/del;
			mcopy(ppart[pp_idx], ppt);
		}
	}else if(m == 1){ /* L Part */
		/* Calculate the L value and put to ppt */
		aurif(b, n, m, ppt);
		for(i=0; i<dcta; i++){
			t = abs(dta[i]);
			if(t == 0)
				continue;
			pp_idx = (t - tbl_first)/del;
			if(dta[i] < 0) {
				// dta[i] has M
				mcopy(tmp, ppartM[pp_idx]);
			} else {
				// dta[i] has L
				mcopy(tmp, ppartL[pp_idx]);
			}

			div3(ppt, tmp, quot, rem);
			if(compare2(rem, 0) != 0){
				err_msg = "error dividing cofactor";
				error(err_msg);
			}
			mcopy(ppt, quot);
		}
		if(n <= pplim){
			mcopy(ppartL[(n-tbl_first)/del], ppt);
		}
	}else {	/* M Part */
		/* Calculate the M value and put to ppt */
		aurif(b, n, m, ppt);

		for(i=0; i<dcta; i++){

			t = abs(dta[i]);

			if(t == 0)
				continue;

			pp_idx = (t - tbl_first)/del;
			if(dta[i] > 0) {
				// dta[i] has M
				mcopy(tmp, ppartM[pp_idx]);
			} else {
				// dta[i] has L
				mcopy(tmp, ppartL[pp_idx]);
			}

			int p;
			set3(&p,tmp);
			div3(ppt, tmp, quot, rem);
			if(compare2(rem, 0) != 0){
				err_msg = "error dividing cofactor";
				error(err_msg);
			}
			mcopy(ppt, quot);
		}
		if(n <= pplim){
			mcopy(ppartM[(n-tbl_first)/del], ppt);
		}
	}
	/* check for intrinsic factor (at most one prime)*/
	t = gcd1(div1(ppt, n, tmp), n); // t = gcd(ppt,n)
	if(t > 1){	/* if we have an intrinsic factor */
		col = col + idec(t)+ 2;
		putc('*', stdout);
		putc('.', stdout);
		j = div1(ppt, t, ppt);	/* divide t out of ppt */

	}
	if(n < start){	/* if n is too small to bother */
		if(n < 40)
			putdec(STDOUT, ppt, 2);
		else{
			putc('o', STDOUT);
			putc('m', STDOUT);
			putc('i', STDOUT);
			putc('t', STDOUT);
		}
	} else{ /* Do some trial division and output to appc */
		if(n ==1){
			while( mod2(ppt, 2) ==0 ){
				col = col+idec(2)+1;
				putc('.', STDOUT);
				div2(ppt, 2, ppt);
			}

		}
		k = n+n;

		for(t=k+1; t<tdlim; t=t+k){
			while(div1(ppt, t, tmp) == 0){
				col = col +idec(t)+1;
				putc('.', STDOUT);
				mcopy(ppt, tmp);
			}
			if(compare2(ppt, 1) == 0) break;
		}
		for(i=col+1; i<74; i++)
			putc(' ', STDOUT);
		putc('C', STDOUT);
		/* k is size of number after parentheses */
		k= mpitoc(ppt, &bigdig, MPCDIM*10);
		i= idec(k);
		apcdec(k);
		apcdec(b);
		fputc(',',file);
		apcdec(n);

		if(m==0){
			if(pm==1)
				putc('+', file);
			else
				putc('-', file);
		}else{
			if(m == 1)
				putc('L', file);
			else
				putc('M', file);
		}
		putc(' ', file);
		putc(' ', file);

		for(j=0; j<k; j++){
			if(j>1 && (j%100 == 0)){
				putc(NEWLINE, file);
				apcdec(k);
				for(i=0;i<12;i++){
					putc(' ', file);
				}

			}
			putc(bigdig[j], file);
		}
		putc(NEWLINE, file);
		fflush(file);
	}
	putc(NEWLINE, STDOUT);
	clear(ppt);
	clear(tmp);
	clear(quot);
	clear(rem);

}

// Find formula for case 2
void aurif(int b, int n, int m, INTEGER ppt) {

	int h, k;
	INTEGER w, x, y, z;
	init(w);
	init(x);
	init(y);
	init(z);
	char *err_msg;

	h = n/b;
	k = (h+1)/2;

	switch(b){

		case 2:
			/* aurif for b=2; n=2h and h=2k-1 */
			mpow(2, h, x);
			mpow(2, k, y);
			add2(x, 1, x);

			if(m == 1)
				sub(x, y, ppt);
			else
				add(x, y, ppt);
			break;

		case 3:
			mpow(3, h, x);
			mpow(3, k, y);
			add2(x, 1, x);

			if(m == 1)
				sub(x, y, ppt);
			else
				add(x, y, ppt);
			break;
		case 5:
		case 6:
			mpow(b, h, z);
			mcopy(x, z);
			add2(x, 1, x);
			mul(x, x, y);
			add(y, z, y);
			mpow(b, k, z);
			mul( x, z, w);

			if(m == 1)
				sub(y, w, ppt);
			else
				add(y, w, ppt);
			break;
		case 7:
			mpow(7, h, w);
			add2(w, 1, w);
			mul(w, w, x);
			mul(w, x, ppt);
			sub2(w, 1, w);
			sub(x, w, x);
			mpow(7, k, y);
			mul(y, x, z);
			if(m == 1)
				sub(ppt, z, ppt);
			else
				add(ppt, z, ppt);
			break;

		case 10:
			mpow(10, h, z);
			mul2(z, 5, x);
			mul(z, z, y);
			mul(y, y, w);
			add(y, z, z);
			add(z, z, z);
			mul2(y, 7, y);
			add(x, y, y);
			add(w, y, y);
			add2(z, 1, z);
			add2(y, 1, y);
			mpow(10, h*3, w);
			add(w, z, z);
			mul2(w, 5, w);
			add(w, y, y);
			mpow(10, k, w);
			mul(w, z, x);
			if(m == 1)
				sub(y, x, ppt);
			else
				add(y, x, ppt);
			break;
		case 11:
			mpow(11, h, z);
			mul(z, z, y);
			mul(y, y, x);
			add(x, z, x);
			mul2(x, 5, z);
			sub(x, y, x);
			sub(z, y, z);
			add2(x, 1, x);
			add2(z, 1, z);
			mpow(11, 3*h, w);
			sub(z, w, z);
			add(x, w, x);
			mpow(11, 5*h, w);
			add(z, w, z);
			mpow(11, k, w);
			mul(w, x, y);
			if(m == 1)
				sub(z, y, ppt);
			else
				add(z, y, ppt);
			break;
		case 12:
			h = n/3;
			k = (h+2)/2;
			mpow(2, h, x);
			mpow(3, k, y);
			mul(x, y, z);
			mpow(12, h, y);
			add2(y, 1, y);
			if(m == 1)
				sub(y, z, ppt);
			else
				add(y, z, ppt);
			break;
		default:
			err_msg = "aurif not defined for b = "+ b;
			error(err_msg);
	}
	clear(w);
	clear(x);
	clear(y);
	clear(z);
}
/* Output n to STDOUT(main) and count columns used */
int idec(int m) {

	char *c;
	int i, nd;
	nd = itoc2(m,&c,10,79);
	for(i=0;i<nd;i++) {
		putc(c[i], STDOUT);
	}
	return nd;
}

/* output k to appc in i4 format
 * k = length of num after ')' */
void apcdec(int k) {

	int i, j;
	char *lildig;
	i = itoc2(k, &lildig, 10, 10);
	for(j=i; j <= 3; j++){
		putc(' ', file);
	}
	for(j=0; j<i; j++){
		putc(lildig[j], file);
	}
}


void initialize2(int t) {
	b = data[t][0];
	n = data[t][1];
	pm = data[t][2];
	del = data[t][3];
	LM = data[t][4];
	LMmod = data[t][5];
	LMres = data[t][6];
	start = data[t][7];
	limit = data[t][8];
	LMlim = data[t][9];

	tdlim = 100000;
	pplim = 1 + LMlim/3;
	tbl_first = n;
	won = pm;
	printf("\nFactorizations of %d^n ",b);
	if(pm >= 0)
		printf("+ %d, n <= %d\n",pm, limit);
	else
		printf("- %d, n <= %d\n", abs(pm), limit);
}

void initialize() {
	// Read xcun.d file
	char *line;	// length of a line in xcun.d
	line = (char *)malloc(sizeof(char)*200);
	int rs = getlin(STDIN,&line);
	int d = 0, i = 0, j=0;
	char *pch;
	while(rs != ERROR && rs != EOF) {
		pch = strtok (line,",");
		while (pch != NULL) {
			d = atoi(pch);
			data[i][j++] = d;
			pch = strtok (NULL, ",");
		}
		i++;
		j = 0;
		rs = getlin(STDIN,&line);
	}
	ntimes = i;

	file = fopen(APPC_FILE,"w+");
	for(i=0;i<1000;i++) {
		init(ppart[i]);
		init(ppartL[i]);
		init(ppartM[i]);
	}
}
